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Abstract 



I The time-dependent form of Tappert's range refraction parabolic equa- 

tion is derived and proposed as an artificial boundary condition for the wave 
equation in a waveguide. The numerical comparison with Higdon's absorb- 
ing boundary conditions shows sufficiently good quality of the new boundary 
^> ' condition at low computational cost. 

1 Introduction 

od _ 

O I In this paper we consider wave propagation in a two-dimensional waveguide fl = 
^ I {{x,y) \ a < y < b} described by the wave equation 

■ T -^^_^_^-n (T\ 

where the sound speed c is a function of spatial variables, c = c{x, y) . For various 
purposes one needs to obtain one way or unidirectional wave equations, which permit 
wave propagation only in certain directions along x-axis. For example, for such 
equations we can consider x as the evolution variable and set initial boundary value 
problems with initial data at x = Xq, which are not well posed for the wave equation 
itself. 
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The main approach to this problem consists in factorization of the wave operator 
in eq. ([T]) into unidirectional pseudo differential factors 

L^L^.L-, L^-^±D. Z, ^ ^(-L^ - |,) . (2) 

where the positive square root is used. Then various approximations to the square 
root are used for obtaining concrete unidirectional equations. Note that the fac- 
torization in eq. ([2]) is not exact when c depends on x, so this approach requires 
optional assumptions of the asymptotic character on this dependence. In fact, 

L = -L- + [D, d/dx] =L+ -L' -D^, 

so it is sufficient to assume that c depends on x slowly, c = c{6x, y) for some small 
parameter 5. Then = 0{6). In this paper we also assume that c is regular 
enough with respect to y to perform the calculations below. 

In the approximation of the square root in eq. ([2]) the assumption of small 
angle propagation with respect to the x-axis is usually used, i. e. the operator 
(l/c^){d'^ /dt^) = A is treated as 0(1) and {d'^/dy'^) = eB is 0(e) in some small 
parameter e. In the case of constant sound speed c, when the operators A and 
B commute, the application of usual rational approximations to the square root 
function then yields local (differential) unidirectional wave equations. Under the 
assumption that X = A + d'^/dy'^ — K is small for some constant the usual 
rational approximations with respect to X can also be applied to approximate the 
square root operator and its exponential (one-step propagator). This approach is 
well developed for the time-harmonic case [U El |3]. 

With the strong dependence of c on ?/ the operators A and B do not commute 
and the applicability of the standard rational approximations to the square root 
operator in eq. ^ is unclear. In fact, difficulties arise already in consideration 
of the first order (linear in B) approximation. As we know now, this problem in a 
particular case was first solved correctly by R. Feynman, who in his paper [1] derived 
the formula ^ 

Jo 

In the late 1970s, F. Tappert in his paper [5j derived the so called range refraction 
parabolic equation for the time-harmonic sound propagation in waveguides with 
arbitrary dependence of the index of refraction on depth: 

^+7)2 f^- ^11^ = 0, (4) 

where n is the index of refraction and ko is the reference wave number. He used the 
expansion of the operator square root in the form very similar to eq. (see eq. ([7j) 
below), but without any derivation. As Feynman's paper was also not cited, the 
derivation of eq. (jl]) in Tappert's paper seemed to be somewhat mysterious. 
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In the next section we derive the nonstationary form of the Tappert's equation, 
using the same expansion, but derived from some elementary results of noncommu- 
tative analysis [H]. We hope that this information will be useful for the reader. 

Note that Tappert's equation (jlj) is of restricted use for the computation of the 
wave field because it is not an amplitude equation. The same is true for our new 
eq. (11 01) . As an application of the new equation we consider its use as an artificial 
(absorbing, non- reflecting) boundary condition (see, e.g., |H1 El El [121 [13]) at the 
boundaries x = Xq. These artificial boundaries and the corresponding boundary 
conditions are needed for restricting the computational domain to simulate waves 
propagating freely in the whole waveguide. Such an use is possible because any 
unidirectional equation partially annihilates waves propagating in the opposite di- 
rection. 

Note that an analogous problem where the variability of coefficients is also es- 
sential was studied by a different method in [HI |T3] . 

2 Derivation of the time-dependent form of Tap- 
pert's range refraction parabolic equation 

Due to the equality 

/ f {sx+{l- s)y) ds = 

Jo x-y 

the first order term in eq. ([3]), denote it C, can be written in the form 

3 1 

2 e^-e^ 

C = eB-^ T' 

A- A 

where the numbers above operators (now called the Feynman numbers) define the 

2 112 

order in which they are to operate [6]. "Thus, BA may be written BA or Ai?" 
jl] . The same formula holds for some general class of operator functions / 

1 3 

f{A + eB) = f{A) + e B ^^^] ~ {^^^ + 0{e^) , (5) 

A- A 

as was proved few years later by Daletskiy and Krein (see [B]). 
Using eq. ([5]), the operator C in the expansion 

{A + eBf" = A^'^ + eC + 0{e^) (6) 
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can be directly computed: 



1 1/2 3 1/2 
2 A — A 2 

C =B ^ =B 



1 3 1 1/2 o 1/2 

^-^ A +1 (7) 

2 /-OO 1 1/2 3 1/2 fOO 

--B / e-^^ e-^^ ds= / e-^^'''Be-^^'''rfs. 



This formula can be verified by squaring both sides of eq. ([6]), which gives the 
equation for C 

A'/^C + CA'/^ = B. 

Substitution into the last formula the expression for C from eq. (JTj) immediately 
gives 

A'/'C + CA'/' = - ^ (e-'^'^' Be-'^'^') ds = B . 

The last equality and the representation (JTj) holds because 

-s/ii/2 f I d\ , 1 . 

e M = exp I 1 u = n(t - , (8) 

and the functions under consideration vanish with their derivatives before some time 
moment, say, for t < 0. 

By a straightforward calculation using eq. ([H]) we first obtain 

r-OO \2 



Cu= / ^.s'-Uuit-2-s,x,y)ds 
Jo c c 







+ / ■s-Ut{t-2-s,x,y)ds 



+ / 'Uyy{t — 2-5, X, y) ds 



(9) 







1 c 
+ 2/ Uyt{t - 2-s,x,y) ■ ds . 
c 



Then, using the formula 

1 c d 1 

- 2-s, X, = -7^^«(^ - 2-s, X, 
c 2 OS c 

and the analogous formula for Uu-, we eliminate the time derivatives and s in eq. 
by integration by parts as follows 

■ utt{t - 2-s,x,y) ds 



r-oo ^ 

c / s ■ Ut(t — 2-s,x,y) ds 
Jo c 

— / u{t-2-s,x,y)ds. 

^ .In C 



At last, changing the variable s by Snew = '^Soid/c, we find the final expression for 
the operator C and so the time-dependent form of Tappert's equation for right- and 
left-propagating waves 



1 If cl\ r , 

u^±-UtT^\Cyy- j / u{t-s,x,y)ds 
1 f°° 

^2 J icuy{t - s, X, y))y ds = 0. 



(10) 



Differentiating eq. (JTOll with respect to t using the above technique, we obtain 
its differential form 

1 1 / 1 
u^t ± -utt T \ Cyy - 7 ) « T 2 = , (11) 

which is, of course, not completely equivalent to eq. ffTUj) . especially in computational 
aspects. Eq. ( fTOl) is more easily handled and, due to its nonlocahty, has better 
stability properties than eq. ( fTTj) . 

It is easy to check that for the time harmonic u = e"^*f(x, y) eq. ( fTTl) coincides 
with the original range refraction parabolic equation for the field v Conversely, 
the Fourier transform with respect to time of eq. (jl]), multiplied by e"^*, gives eq. fill I) 
and then, after integration with respect to t, eq. (fTOl) . 



3 Eq. (1 101) as an artificial boundary condition 



In this paper we will apply eq. fllOj) for left going waves as an artificial boundary 
condition on the artificial boundary x = xq, the computational domain is located in 
X > xq- As was mentioned in the Introduction, the problem consists in simulating 
waves, propagated in the unbounded waveguide, by solving some initial boundary 
value problem in the restricted computational domain. 

As far as we know, only the technique of perfectly matched absorbing layers can 
be immediately applied to the problem under consideration in the case of waveguides 
with variable sound speed c (see [8]). Besides this technique the most frequently 
used by the practitioners boundary conditions have the form 



n 



dt ^ dx 



u = 0, (12) 



where Cj are some constants. In the paper Higdon proved (in the case of con- 
stant c) that the class of these conditions contains all boundary conditions obtained 
by the Fade approximations of the operators L"^ (the Engquist-Majda boundary 
conditions [7]). It is immediately seen that the condition (fT2l) perfectly annihilates 
waves of the form f{Ct — x)(l){y) if the phase speed C is equal to some Cj. The 
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boundary conditions (fT2|) are referred in the literature as the Higdon conditions of 
order J. 

The main problem in any practical use of the Higdon conditions consists in 
finding some (quasi) optimal set {Cj}. For waves with given frequency spectrum it 
is possible to try the phase or group velocities of the corresponding normal modes. 
The automatic choice, described in [19] , also gives good results for dispersive waves. 
In essentially non- stationary problems for the wave equation the correct choice of 
Cj is difficult and the optimal set {Cj} may depend on time. So in such ambiguous 
situations the use of some mean value of the sound speed c as Cj is probably more 
preferable. 

As our boundary condition reduces to the second order Engquist-Majda condition 
in the constant sound speed case, we expect that it will work better than the Higdon 
condition of order 2 in the variable sound speed case. The performed numerical 
experiments confirm this conclusion. 

We must make some remarks on the well-posedness of the introduced boundary 
condition. In the case of constant sound speed one can recognize in eq. (11 II) the 
second order Engquist-Majda boundary condition [7], which is proved to lead to a 
well-posed mixed initial boundary value problem for the wave equation. The proof 
consists in checking the uniform Kreiss condition (U.K.C.) and is valid as well in the 
variable coefficient case (see, e.g., ^J). As the third term in eq. ( ITTl) is unessential 
for the U.K.C, we may conclude that the boundary condition ( ITTl) leads to a well- 
posed mixed initial boundary value problem for the wave equation. As eq. (ITOl) 
implies eq. (fTTI) . the same is true for the boundary condition ( ITUl) . 

4 Finite difference discretization 

In our numerical experiments we have used the standard second-order explicit finite 
difference scheme for the wave equation on the uniform grid {t\Xj,yi} 

^D+D;uli = DtD-uli + D^D-uli, (13) 

where i = u{f, Xj, yi){f' = ir, Xj = jh — h/2, yi = lh,i = 0,1 j = 1,2 ... ,1 = 
1,2...), Dtu], = (4+1^, - ul,)/h, D~ul, = (n}_, - n}_,,,)//i, and D+, Dr, , 
are defined analogously. 

We assume that u and all its derivatives vanish for alH < 0. Under this assump- 
tion the boundary condition of the Tappert type (fTUj) on the left boundary at x = 
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Figure 1: Computational domain for the point source. 



are discretized as: 



DUul, + </)/2 - Dtiul, + «^,,)/(2c3A,) + 

° \ ^ J 1,12.1 ^-9 



3/2,j rn=2 
i 

m=2 



(14) 



m=2 



where (■)3/2j is the value of (■) at a; = 0, y = jl. It should be noted that in practice 
the sound speed c is often given by some interpolation formula, so there is no need to 
calculate its derivatives by the difference methods. The Higdon boundary conditions 
were discretized as described in [18]. The implementation of the boundary condition 
(imp in the finite element framework can also be easily obtained. 



5 Numerical experiments 

The nondimensional computational domain for these experiments is shown in fig. 1. 
It is a rectangle {x,y\ < x < = 10 ; < y < Ly = 10}. At each boundary of 
the extended domain A'BCD' the hard wall (zero Neumann) boundary conditions 
are applied but waves do not reach these boundaries during calculations. At the 
boundary AD of the truncated domain ABCD the absorbing boundary conditions 
are imposed. Thus solving the initial boundary value problem for the wave equation 
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Ly/2 



Ly 



Ly/5- 



Figure 2: c(y) given by (|T71) . 



Figure 3: c(y) given by (|T8 



first on the truncated domain and then on the extended one we can estimate the 
quahty of solutions using error measures 



'EJ2\iJ'ik-UE{ti,Xj,yk)\ 



and 



e{ti 



Y.Y.\uE{ti,Xj,yk)\ 

j k 



max max \u) — UEiti, xj, yk)\ 

j k 



(15) 



(16) 



where u is the solution of initial boundary value problem on the truncated do- 
main, and M^; is the solution from the extended domain. 

As the initial conditions was taken the wave field in the homogeneous media 
produced by the point source of short duration located at the center of the extended 
domain. 

Numerical solutions were obtained on the grid with the space step a = 0.1 and 
the sufficiently small time step r which guarantees stability of the finite-difference 
scheme. 

In our experiments we compare the new boundary condition, in the sequel called 
the Tappert condition, with the Higdon boundary conditions of order J = 2 and 
J = 3, with the {Cj} given by 



^q" c{y) dy 



For the first example we use c(|/) with a minimum at the depth Ly/2 given by 
(fig- 2) 



CUE 



0.5e^ 



(17) 

which gives Cj ~ 0.8565. The point source duration was taken to be 1.4. The 
results of calculations are presented in fig. 4-5. In the second example the variation 
of c(?/) is larger (fig. 3), where ciy) is defined through 
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^ Tapped 

■ - Higdon, order 3 

Q- Higdon, order 2 




5 2 2.5 3 3.5 4 4.5 5 5.5 




Figure 4: Errors E{t) for the first ex- Figure 5: Errors e(t) for the first ex- 

ample ample 

y 

c{y)=4-^ [ e^'-'^y^'^' ds . (18) 
V'^ J 

— oo 

In this example Cj ~ 1.6521 and the point source duration was taken to be 1. The 
corresponding results are presented in fig. 6-7. We see that the larger the variation 
of the sound speed is, the better the Tappert condition works. 

As the Tappert condition ffTUj) contains no derivatives of the sound speed with 
respect to the range variable, it can be applied in the range-dependent environment 
without any modifications. If the third example the following range- dependent sound 
speed c{x, y) was taken: 

c(x,?/) = l-0.5e("-°■^^-)'/^ (19) 

In this example the source was positioned at a; = 3Lj./4, y = Ly/2. The results 
are presented in fig. 8-9. This and various other experiments show that in the range- 
dependent case the growth of errors under the Tappert condition ( ITOi) is essentially 
the same as in the range-independent examples. 

The discrete form of the Tappert boundary condition 0141) shows that the inte- 
grals in the eq. (fTUj) are simply accumulated at each time step, so the computational 
cost of the new boundary condition only slightly differs from those of the first or- 
der Higdon condition. The computational cost of the second order and the third 
order Higdon conditions is about 15% and 40% greater. These estimates are very 
approximate because the calculations were conducted in the MATLAB framework. 

At last we should note that the behavior of the new boundary condition over 
long times is essentially as good (or bad) as those of the Higdon conditions. 
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^— Tapper! 

- - Higdon, order 3 

-0- Higdon, order 2 



10 ' - 




Figure 6: Errors E{t) for the second 
example 



^— Tappert 

— - Higdon, order 3 

-e- Higdon, order 2 



10"^ r 




0.5 1 1.5 2 2.5 3 



Figure 7: Errors e(t) for the second 
example 



6 Conclusion 

In this paper the time-dependent form of Tappert 's range refraction parabolic equa- 
tion ffTOl) is derived. This equation describes unidirectional, small angle (with respect 
to the X-axis) wave propagation and can be used for solving various problems of non- 
stationary wave propagation which are, in particular, ill-posed for the wave equation 
([T]). Despite its nonlocality in time, it is easily calculated. 

Here the new equation was applied as an artificial boundary condition for the 
wave equation in a stratified waveguide. Numerical experiments, performed for 
comparison of the new boundary condition with Higdon's absorbing boundary con- 
ditions, showed its sufficiently good quality at low computational cost. 

Eq. (ITOl) . regarded as an artificial boundary condition, can be considered as the 
first member of a new sequence of boundary conditions which are suitable for strat- 
ified waveguides. This sequence can be obtained by the calculation of the following 
terms in expansion of the square root operator into the Newton series (noncommu- 
tative analog of the Tailor series, see [H]). To achieve good stability properties the 
regularization of obtained in such a way partial sums will be necessary. This can 
be done by the application of some form of noncommutative rational approximation 
(e. g.. Fade appoximation [20l Fart 2, Sec. 1.5]) and will lead also to the wide angle 
extensions of eq. (ITOl) . 

The problems for further investigations consist in the calculation of the following 
members of this sequence, verification their well-posedness and automatization of 
their numerical implementation. We intend to consider some of these problems in 
our future studies. The results of our future studies will be compared with the 
results of the recent works flU{ [TT] on the Higdon boundary conditions, where the 
problem of variable coefficients was also considered. 
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xio"^ Sum-error E(t) ^^g-s Max-error e{t) 




Figure 8: Errors E[t) for the third ex- Figure 9: Errors e{t) for the third ex- 

ample ample 
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